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ABSTRACT 

This  paper  studies  prediction  of  future  failure  (rates)  by  hierarchical  empirical  Bayes  (EB)  Poisson 
regression  methodologies.  Both  a  gamma  distributed  superpopulation  as  well  as  a  more  robust 
(long-tailed)  log  student-t  superpopulation  are  considered.  Simulation  results  are  reported 
concerning  predicted  Poisson  rates.  The  results  tentatively  suggest  that  a  hierarchical  model  with 
gamma  supcrpopulation  can  effectively  adapt  to  data  coming  from  a  log-Student-t  superpopulation 
particularly  if  the  additional  computation  involved  with  estimation  for  the  log-Student-t 
hierarchical  model  is  burdensome. 


1.  INTRODUCTION 

The  following  model  often  provides  a  useful  place  from  which  to 
commence  the  analysis  of  point  event  process  data.  First,  suppose  there  is  a 
set  of  I  entities  or  units,  each  of  which  generates  an  observed  history  of  point 
events.  Take  each  describing  point  process  to  be  homogeneous  Poisson  (\j),  i 
=  1,  2,  3,  ...,  I.  The  observed  data  appears  as  (s;,tj),  s*  being  the  number  of 
events  for  process  i  over  active  or  operating  time  tj.  Also  observed  are  certain 
fixed  explanatory  variable  values;  x*j,  j  =  1,  2,  ...,  p,  associated  with  A,j.  In  some 
literature,  e.g.,  Everitt  (1984),  such  variables  are  called  manifest.  Second, 
there  is  a  latent  quantity,  5;,  associated  with  that  is  unobservable  but 
influences  X\  behavior.  It  is  convenient  to  view  8j,  at  least  provisionally,  as 
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bein3  drawn  randomly  from  some  superpopulation  of  values  and  held  fixed 
thereafter,  thus  endowing  Xj  with  its  own  particular  individuality. 

We  call  such  a  setup  hierarchical,  and  ask  it  to  furnish  insights  and 
numbers  concerning  (a)  the  individual  rate  values,  Xj,  (b)  the  influence  of  the 
explanatory  variables  upon  these  rates,  and  (c)  the  nature  of  the 
superpopulation  that  gives  rise  to  the  latent  variable  values;  future  values  of 
the  rates,  e.g.,  Xi+i,  etc.,  may  be  viewed  as  coming  from  such  a  population,  at 
least  to  a  first  approximation. 

The  above  model  suggests  itself  for  many  purposes,  one  in  particular 
being  in  risk  analysis,  e.g.,  of  nuclear  power  plant  safety  systems.  Such  setups 
are  also  natural  in  other  reliability-related  areas  as  well,  particularly  in  ones 
arising  in  the  military.  Application  may  perhaps  be  made  to  data  reflecting 
"human  unreliability,"  i.e.,  the  propensity  of  different  individuals  to  make 
errors,  or  experience  accidents. 

The  purpose  of  this  paper  is  to  describe  methods  for  fitting  various 
hierarchical  models  to  the  type  of  data  described.  Particular  attention  is 
devoted  to  the  prediction  problem:  given  the  past  record  of  an  individual 
item  (e.g.,  human  being),  how  well  can  one  predict  its  (her)  future 
performance,  even  if  some  basic  conditions  change? 

2.  THE  FORMAL  MODEL 

The  formulation  proposed  can  be  written  as  follows:  for  i  =  1,  2,  ...,  I,  and 

fl=(Pl,P2,--,Pp)T 

Sj~!ID  g(-;fi) 

Xj  =  f(xip,6i),  (2.1) 

Sj  I  Xj,tj~Ind.  Poiss(Xjtj); 

fi  =  (0i,  02,  -,  ©r)  being  a  parameter  identifying  g,  the  density  associated  with 
the  assumed  fixed  superpopulation.  In  what  follows  we  concentrate  on 
certain  parametric  forms  for  the  link  function  f  and  the  superpopulation  g, 
and  aim  at  estimating  the  ft-value  best  representing  the  superpopulation 
giving  rise  to  the  apparent  X-values.  For  various  reasons,  convenience  and 
tradition  being  influential,  we  restrict  attention  to  the  log-linear  model 

X,  =  f(xj£,5j)  =  exp(xjj}+5j).  (2.2) 

As  suggested  earlier,  the  objectives  of  the  analysis  will  be  several-fold,  but  an 
important  one  will  be  to  estimate  an  individual  Xj-value,  i.e.,  the  actual 
realization  of  Xj  that  prevails.  An  even  more  important  objective  is  to  predict 


2 


the  number  of  future  events  associated  with  i,  Sj(t).  This  entails  finding  an 

A  A 

estimate  Jl,  and  one  for  the  individualizing  parameter  5j,  namely  8;. 
Estimation  will  be  carried  out  by  assuming  that  g(-,Q),  the  superpopulation 
density  giving  rise  to  5*,  is  one  of  a  specific  parametric  family,  and  first 
estimating  the  parameters  of  that  density  along  with  the  regression 

A  A 

parameters.  At  a  later  stage,  the  estimated  parameters  Ji  and  &  are  utilized  to 
create  estimates  of  8j,  and  finally  X;,  see  Cox  and  Hinkley  (1974),  p.  401,  Morris 
(1983),  Deely  and  Lindley  (1981),  etc.  The  several-stage  or  hierarchical  analysis 
is  referred  to  as  parametric  empirical  Bayes  (PEB). 

This  work  is  an  extension  of  Gaver  and  O'Muircheartaigh  (1987)  in  which 
discrepancy-tolerant  (robust)  estimates  of  8j  and  X*  were  produced  and 
evaluated  without  consideration  of  explanatory  variables.  The  major 
purpose  of  the  present  article  is  to  consider  the  effect  of  explanatory  variables 
in  the  context  of  hierarchical  models  using  quite  different  models  for 
superpopulations:  first,  the  simple  conjugate  Gamma,  and  next  the  log- 
Student  t  with  a  small  number  of  degrees  of  freedom  so  that  tails  are 
extended,  and  outliers  more  apt  to  be  generated. 


3.  AN  EMPIRICAL  BAYES  APPROACH 

The  approach  taken  to  providing  estimates  is  traditional;  see  Berger  (1985, 
Chap.  4).  We  first  remove  the  condition  on  5  for  each  item  to  obtain  the 
unconditional  likelihood 


L(^,£;s,f)  =  nJe 

1  =  1 


■'fet'fc  (/(*,&*)>, f 


S.! 


g(5;i)dS. 


(3.1) 


A  A 

The  latter  is  then  maximized  with  respect  to  £  and  Q.  to  produce  JJ  and  ft. 
These  quantities  are  then  inserted  into  the  expression  for  the  posterior 
density  of  8, 


gr(S)~  gr(S;e,s„t„  xj)  =  Kle-'kl,>''(f(xle,s)l,f,g(S;e) 

where  the  constant  Kj  is  a  normalizing  factor.  A  point  estimator  of  Xj  is  taken 
to  be  the  posterior  mean  (other  options  of  course  available). 


(3.2) 
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where  x,  is  the  value  of  the  explanatory  variable  for  conditions  anticipated 

when  the  estimate  is  to  be  applied;  if  x,  =  x,  then  we  have  Xit  an  empirical 
Bayes  estimate  of  A.},  for  conditions  under  which  the  data  were  taken;  this  will 
often  be  a  shrunken  estimator  that  has  a  smaller  mean-squared  error  than 
does  a  simple  individual  estimator.  If  x,  refers  to  other  (e.g.,  future) 

conditions  then  A,  calculated  by  (3.2)  may  be  called  the  mean  predictive  rate. 
If  more  information  is  desired  then  the  entire  predictive  distribution  is 
needed: 


m)  =  J e',M  *  h'W* 


(3.3) 


this  approximates  the  conditional  probability  of  s{  future  events  for  item  i, 

given  that  it  is  exposed  for  time  f,  and  under  conditions  x,. 

It  is  apparent  that  the  approximation  so  obtained  may  be  under-variable, 

in  that  it  treats  /3  and  0  as  fixed  and  known  in  (3.2)  and  (3.3).  The 

hierarchical  Bayes  analysis  described  by  Berger  (1985,  Chap.  4)  is  a  substitute 
that  avoids  that  criticism.  This  defect  is  undeniable,  but  some  appreciation 
for  the  magnitude  of  the  effect  can  be  obtained  by  bootstrapping.  Of  more 
concern  to  us  has  been  investigation  of  the  effect  of  superpopulation  model 
choice:  how  different  can  actual  prediction  be  in  simple  situations  modeled 
quite  differently?  We  proceed  to  compare  and  contrast  two  models,  one 
conjugate  Gamma  and  the  other  longer-tailed  and  hence  outlier-prone. 


4.  GAMMA  LATENT  VARIABLE  POISSON  REGRESSION  (GALVPR). 

It  is  conventional  and  convenient  to  invoke  the  gamma  density  to 
represent  the  random  effect  in  (2.1);  see  Lawless  (1987  a,b)  and  Anscombe 
(1950)  for  examples.  Thus 


g(M 


- u/a(u/af 1  1  1 

rf"  ')  ‘ 


is  the  superpopulation  model,  from  which 

E[S]  =  1 
Var[8]  =  a. 


(4.1) 


(4.2) 
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Lawless  (1987b)  gives  expressions  for  the  ln-likelihood  and  its  derivatives  for 
this  hierarchical  model.  It  turns  out,  however,  that  a  more  satisfying 
parameterization  is  in  terms  of  0=ln  a  when  the  mle  stage  is  undertaken. 
Since  a  one-parameter  gamma  density  is  used,  the  regression  has  a  constant 
term;  that  is  xjj  =1.  For  convenience  we  provide  expressions  for  the  ln- 
likelihood  and  its  derivatives  using  our  parameterization. 

In  the  present  parameterization,  then, 

X  =  U  exp[x$]  (4.3) 

so  exp(S)  =  U  is  gamma.  In  order  to  form  the  likelihood  element  in  (3.1)  it  is 
only  necessary  to  integrate  to  obtain  the  explicit  form 


ila  H  e  (Unltfili 

Us,r(£- 


Xs5a_rf_!_y 


where  c,  -  expjx,/}}.  The  log-likelihood  is 


L 


/(e,^;s,f)  =  X  +  -  In  Si !+  s,|ln  qtj  -  ln|edc,tj  +  l)j  -  e  e]n{eeciti  +  \)j 

,=i  u=° 


where  ^ln(l  +  ie°)  =  1  if  S,  =0. 

,=o 

The  following  derivatives  can  be  obtained. 


/  hi- 1 


Itt1 


i=l  [(  j= 0  1  +  je 


lnle9ci',+1l'7^l[s,+c"el  a6) 
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(4.9) 


d^l  ^  sie  cih+ciU 

[e°cili+lf  v,i 


where  the  summations  involving  Sj-1  are  set  equal  to  zero  when  Si=0. 
Further, 


’  -d7l  ' 

1 

_  v 

«  \ 

dpjp, 

Ldt. 


(4.10) 


A  Newton-like  iterative  procedure  is  used  to  solve  the  system  of 
equations 


^  (Os 

II 

O 

(4.11) 

II 

o 

(4.12) 

If  S{  is  large  then  evaluating  the  sums  appearing  in  (4.5),  (4.6),  (4.7),  (4.13) 
and  elsewhere  tends  to  be  time-consuming.  However  all  such  sums  are  well- 
behaved  (of  monotonic  formations)  and  can  be  well-approximated  by 
integrals.  This  feature  is  not,  but  easily  can  be,  included  in  our  programs. 

If  (3k)  were  known,  then  a  Newton  procedure  to  estimate  0  would  be  to 
recursively  solve  the  linear  equation 


dl_ 

99 


f  d2 


e  =  e° 


+ 


9ll 


99 2 


e  =  e° 


(0-9°) 


(4.12) 


n  01 

where  0U  is  a  current  estimate  of  0.  Note  that  if  =  0,  then 


d2l 


de 1 


=  g(0)  =  e 
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/ 

f  s.-l 

(  ■ 

2' 

I 

I- 

i 

,  o 

i=i 

[i=° 

(l  +  /e  ) 

•  2e 


-39,  r  0, 


IneV.+l  +-J 


2cLt,e 


-20 


eV,+l 


c,‘, 


e  c,t,  + 1 ) 


[s,  +e  e] 

(4.13) 


Hence,  (4.12)  can  be  rewritten  as 


6 


e-eu  = 
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i= 1 


+e_20HA^+l]--e£^—[s,  +  ee]- 
;=0  ^  +  ft  e  ci^i  +  ^ 


/ 

f  *1-1 

(  i  '\ 
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1  •  6 
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[\  +  je  J 

/ 
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-26  ( 


ec,t,  +1 
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(4.14) 


To  obtain  an  initial  estimate  of  0,  note  that,  letting  Nj(tj)  denote  the  ith 
random  variable  of  the  number  of  observed  events, 

E[N,{t,)]  =  C,t,;  (4.15) 

Var\Nl{t,)]  =  c,tl[U-c,tlee}.  (4.16) 

Thus,  [(Ni(ti)  -  Cit^/Vcib]  has  mean  0  and  variance  [1  +  Citje9].  We  propose 
starting  the  iterative  procedure  to  find  0  by  computing 


(4.17) 


If  m  <  then  a  log-linear  model  is  used  to  describe  the  data.  If  m  >  1,  then 
the  initial  estimate  of  0  is 


So  =ln 


(m-l) 


7**" 

. 1  1=1 


(4.18) 


If  0  were  known,  then  (P)<)  could  be  estimated  with  generalized  linear 
model  software  in  the  following  manner,  (cf.  McCuIIagh  and  Nelder  [ I9S3J). 

A  Newton  iteration  to  solve  the  equations  0  =  ~  is  to  solve  the  system  of 

dpk 

equations 
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dl 
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WJP, 


(0  =  0°) 


(Prt>°)  (419) 


where  is  the  current  estimate  of  {i. 
Put 
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c,t,  7 
e6c,t,  4- 1 


where  q  =  exp{x$°). 

Equation  (4.19)  can  be  rewritten  as 


/  /  p  ^ 

°  =  S  y.-£w.A  u* 

■=iV  >=l  y 


k=l, ...,  p 
where 


y.  =[s,-  c,f,  J[c,i t,(e‘c,t,  + 1)] 


+  I“AC 

7  =  1 
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and 


(4.20) 


(4.21) 


(4.22) 


(4.23) 

The  equations  of  (4.21)  are  the  normal  equations  for  a  least  squares  regression. 
The  following  is  an  iterative  procedure  to  obtain  estimates  of  0  and  (pk) 
a)  Fit  a  log-linear  model  stopping  after  one  iteration 
1.  Start  with 


=  In 


M)/,( 


Z  Solve  the  equation  (4.21) 


with 


iv,  = 


1 

S'  +  2 


(  < 


y.  =  -«>, 


uy 


b)  Find  the  initial  estimate  of  0  by  evaluating  (4.18).  If  m<l,  use  the  log- 
linear  model  of  a)  to  describe  the  data. 

I.  Next  estimate  (pi<.):  Evaluate  and  solve  equations  (4.20)  -  (4  23). 
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II.  Next  estimate  0:  Evaluate  and  solve  equation  (4.14). 

III.  Continue  alternating  between  I  and  II  until  convergence. 

In  the  simulation  experiments  described  in  Section  7  the  above-obtained 
estimate  of  0  occasionally  either  cycled  among  negative  values  or  became 
very  large  and  positive.  In  these  cases  II  was  replaced  by  a  search  of  the 
marginal  likelihood  for  0  with  fixed  (PiJ. 

5.  ROBUST  HIERARCHICAL  POISSON  REGRESSION  (ROLVPR):  THE 

LOG-STUDENT  t  SUPERPOPULATION 

As  an  alternative  to  the  GALVPR  model,  allow  8  to  have  the  Student  t 
density 

g(8;t2,d)  = - (5.1) 

(l  +  <$/f2)  2 

this  distribution  is  adjustably  longer-tailed  than  is  the  log-Gamma 
distribution  (for  6)  of  the  previous  model,  and  hence  better  represents 
outliers  and  extreme  extra-Poisson  variability.  The  parameter  d  is  the 
"degrees  of  freedom"  for  the  Student  t;  for  the  present  purpose  a  low  value  of 
d  (e  g.,  d  =  3-5)  is  useful.  The  Student-t  model  for  log  failure  rate  was 
introduced  in  Gaver  and  O'Muircheartaigh  (1987).  There  it  was  pointed  out 
that  the  marginalization  step  of  (3.1)  could  be  performed  using  Gauss- 
Hermite  numerical  integration;  see  Naylor  and  Smith  (1982).  In  this  paper 
we  employ  a  variant  of  the  Gawss-Hermite  technique  that  involves  an  initial 
correction  by  Laplace's  method. 

The  procedure  currently  adopted  for  fitting  the  regression  parameters  (i 
in  addition  to  the  Student  t  parameter  x  proceeds  iteratively:  first  explain  as 
much  item-to-item  variability  as  possible  by  suitably  weighted  regression, 
then  alter  the  model  to  approximately  adjust  for  regression  effects  and  apply 
the  methodology  of  Gaver  and  O'Muircheartaigh  (1987)  to  estimate  x2.  This 
value  then  provides  refined  weights  for  a  new  regression.  We  speak  of 
rocking  back  and  forth  between  the  regression  and  latent  variable  stages. 

5.1.  Rocking  Algorithm  when  6;  -  Student  (p,  X,  d) 

Here  is  how  the  above  procedure  operates  when  latent  variables  are 
Student  t  so  as  to  represent  adjustably  long-tailed  outlier-prone  regressions;  d 
>  1  is  a  tuning  parameter  with  Var(5]  =  x2d/(d-2)  if  d  >  2. 

a)  Regress  yj(l)  =  Vsi  ln(sj/tj)on  Vs*  xj;  (5.1) 
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Replace  sj/tj  =  0/tj  by  l/3t,.  Obtain  g(l). 

b)  In  the  ith  likelihood  component  obtained  by  integrating  out  with 
respect  to  the  5,-distribution, 

=  ~dz  (5.2) 

[l  +  (z2  /  T 2d)]  2  T 

where 

lnA.(z)  =  x,£  +  Z/  (5.3) 

replace  t;  by  =  f;(l).  Now  numerically  optimize  (5.2)  by  choice  of 

r=f(2);xA(i)  is  a  moment  estimator.  Details  of  the  likelihood  integral 
approximation  and  optimizations  are  furnished  in  Section  6. 


c)  Regress  y.(2)  = 


d  W 


7+?(2).  2 

Vsi 


f  \ 

In  (s/f;)  on  — I-  "P(2) 


d-2) 


where  tj  and  s,  are  the  original  data  values.  Obtain  P(2). 


d)  In  step  b  above  replace  ti  by  f,eiifi<2)  and  x,$  by  x,£(2)  and  again 
numerically  optimize  to  find  t( 3). 

e)  Return  to  step  c)  with  T^p). 


f) 


Continue  to  convergence  of 


{£(*)},  {^(fc)}. 


The  above  procedure  converges  rapidly  in  our  experience,  giving  results 
in  close  agreement  with  the  simultaneous  optimization  of  the  likelihood 
with  respect  to  x  and  p.  The  latter  is  a  much  more  computationally 
demanding  procedure  than  is  rocking. 


6.  LIKELIHOOD  COMPUTATION 

An  essential  part  of  the  preceding  algorithm  is  the  numerical  evaluation 
of  likelihoods  of  this  form: 

l(*1;{s,L))  =  (61) 

i=1 

where 


10 


(6.2) 


[A,(z)r  m  i 


=  J  e~Q,(z)dz. 

Under  the  log-linear  model 

lnA,(z)  =  *J  +  z  (6.3) 


so 

Q,(z)  =  A,  (z)t,.  -s,  In  A,  (z)  +  ^y^-ln[l  +  z2  /  T2d]  +  lnT,  (6.4) 

omitting  irrelevant  constants.  In  order  to  evaluate  the  integral  in  (6.2) 
approximately  but  reasonably  accurately  we  apply  either  (a)  a  version  of 
Laplace's  method,  in  which  Qj(z)  is  approximated  by  a  quadratic  and 
integrated  explicitly;  alternatively  (b)  apply  a  refined  version  of  (a)  involving 
Gauss-Hermite  integration  of  the  error  resulting  from  the  quadratic 
approximation  to  (6.3).  Here  is  a  sketch  of  the  process.  In  what  follows  we 
will  modify  the  time  to  be  t{e^. 

6.1.  Laplace  Method,  and  a  Refinement. 

To  compute  L,  =  j  e~QMdz,  the  ith  likelihood  component,  we  begin  by 

approximating  Q;  by  a  quadratic  as  follows.  Since  A*  =  e2,  t;  is  modified  as 
indicated  above,  and 


-Q,(z)  =  -A.f,  +  s,  In  A,  -  ^Y^ln(l  +  z2  /  r2d)  -  In  r, 
dQ,  (d  +  \\  1  z 

dz  ~  l+S‘  {  i  J[l  +  Z:  /  T2dj  T2 


(6.5) 


=  -elf,  +  Sj  -  w(z)z  /  t2 

where 

w(z)  =  ^-^-[l  +  z2  /  r2d] 
d  1 

is  considered  to  be  a  weight.  Now  -dQj/dz=0  entails  the  equation 
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e 


(6.6) 


r  d+ 1  y 

1 

z 

c  4- _ 

d 

t, 

‘  T2 

[l  +  z2  /  r2d] 

_ 

K  ). 

Equation  (6.6)  may  have  two  solutions.  We  obtain  a  single  reasonable 
approximate  solution  to  (6.6)  as  follows:  An  initial  solution  to  (6.6)  is 


2,(0)  =  ln(s,  /  ft)  if  s,.  >  0, 
z,(0)  =  ln(l/3f,)  if  s,  =  0; 


(6.7) 


other  replacements  for  the  zero  count  situation  are  possible. 

Let  z,  be  the  solution  to  the  equation 

et=j[si-zw(zi{0))/t2] 

after  one  Newton-Raphson  iteration  starting  at  Zj(0). 

Next  evaluate  an  approximation  to  Q,  (z,).  First  approximate 


Qi(2)  -  e\  -  s,  +  (6.8) 

T 

Hence, 

Q[(z)- e**i+U>( 2,)A-  (6.9) 

T 


Finally  approximate  ez‘  by  Sj/tj  resulting  in  the  approximation 

Qi(z,)  ~  si  +w,((2i)^2-  (6.10) 

Write 

Q,(z)  =  Q,(z,)  +  (2  -  z,)2  |q,"(z,)  +  ^,(z);  (6.11) 

Laplace's  method  assumes  that  R,(z)  is  negligible  and  hence 
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(6.12) 


(-.(rVvt.)  =  jyal‘Vz 

-.-'vs/VaW 

so  the  log-likelihood 

t(r2;s,t)  =  “X  [<3,  (2. )  +  ^ ln  Qi "  (z. )  J  (6.13) 

which  can  be  numerically  optimized  by  choice  of  x2  for  fixed  tuning  constant 
value  d  (in  principle  optimization  on  d  can  also  be  included). 

Improved  numerical  results  have  been  achieved  by  writing 

^;a)=r^ 

(6.14) 

=  e-Q(~2-]p  /  Qr(zt)Sy^e-^ 

with 

R*(w)=Ri(-yl2/Qi"(zi)w  +  zi),  (6.15) 

Ri (w)  being  defined  by  (6.11).  The  integration  is  then  performed  by  Gauss- 
Hermite  technique,  i.e.,  by  replacing  the  integral  by  a  finite  sum  at  points  Wj 
determined  by  zeros  of  the  Hermite  polynomials;  see  Abramowitz  and  Stegun 
(1964).  Experience  has  shown  that  the  above  produce  numerical  results  that 
agree  well  with  other  numerical  methods  such  as  that  of  Naylor  and  Smith 
(1982);  however,  the  unadorned  Laplace,  (6.12),  may  sometimes  be  satisfactory, 
and  is  certainly  more  quickly  computed,  which  is  a  virtue  if  bootstrapping  is 
undertaken. 

Alternative  computational  procedures  exist  and  have  virtues.  The 
Newton-like  iteration  applied  to  the  Gamma  model  of  Section  4  can  be 
adapted  to  the  log-Student  model,  but  we  have  not  undertaken  this  as  yet.  A 
sampling-based  approach  of  Gelfand  and  Smith  (1988)  is  a  natural  option,  but 
at  present  appears  unnecessarily  computer-intensive.  As  will  appear,  even 
the  apparently  crude  rocking  approximation  proposed  leads  to  interesting 
contrasts  between  predictions  made  by  the  conventional  conjugate  Gamma 
and  the  robustifying  Student. 
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7.  NUMERICAL  ILLUSTRATIONS 

In  order  to  illustrate  the  performance  of  the  two  proposed  prediction 
schemes  we  have  performed  extensive  simulations.  These  illustrate  the 
anticipated  comparative  performance  of  GALVPR  and  ROLVPR:  the  latter  is 
often  better  able  to  adapt  to  the  appearance  of  large  outlier  rates  by  refusing  to 
shrink  them  down  as  extensively  as  do  the  former.  The  difference  between 
the  predictions  made  by  the  two  schemes  is  less  noticeable  for  small  rates; 
here  the  behavior  of  the  gamma-based  approach,  GALVPR,  may  actually  be 
superior,  probably  because  of  the  approximations  made  when  implementing 
the  Student  ROLVPR  model.  Improvements  in  the  current  procedure  for 
fitting  the  latter,  e.g.,  when  counts  are  zero,  are  likely  to  show  up  as  reduced 
upward  shrinkages. 

Simulation  Experiment 

The  present  simulations  are  all  based  on  a  group  of  1=20  items.  For  the 
log-linear  rate  of  (2.2)  x^  =  pi  +  p2xi/  and  xj  =  +1  for  i  =  1,  2, ...,  10,  Xj  =  -1  for  i  = 
11,  12, ...,  20.  In  addition  pi=0.5  and  (32  =  0.1  and  0.3,  while  tj  =  2  throughout. 
For  each  experiment  20  Poisson  rates  were  then  generated  from  the  Student 
model  with  x  =  1.0  and  d  =  5,  and  for  each  rate  a  single  Poisson  data  point  was 
generated  with  mean  A.,ti.  These  then  constitute  the  observed  counts  from 
which  predictions  are  made.  Each  prediction  is  viewed  as  a  point  estimate  of 
the  underlying  Poisson  mean  giving  rise  to  the  corresponding  observed 
count;  it  is  a  natural  point  estimate  for  a  future  count.  The  predictions  are 
chosen  to  be  the  means  of  the  posterior  distributions  from  the  GALVPR  and 
ROLVPR  model  specifications,  where  each  model  is  fitted  to  the  data  (20 
counts,  plus  values  of  Xj)  for  the  particular  experiment,  meaning  that  (3j,  i  =  1, 
2,  and  a,  for  GALVPR,  and  x2,  ROLVPR  were  estimated  as  described  earlier. 
These  models  were  actually  fitted  by  two  methods:  (a)  to  all  count  data  in  the 
experiment,  including  that  for  the  item  whose  rate  is  predicted,  and  (b)  to  all 
data,  but  omitting  the  observation  for  the  item  to  be  predicted,  i.e.,  in  cross- 
validation  mode. 

An  illustration  of  a  particular  experimental  outcome,  and  the 
corresponding  predictions  appears  in  Table  I.  Note  that  for  this  particular 
data  set  the  average  mean  square  error  of  ROLVPR  no-cross-validation 
predictions  is  the  smallest.  This  is  not  always  so;  see  the  figures  for 
comparisons  of  mean-squared-errors  for  the  two  shrunken  predictions,  and 
raw  predictions. 
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TABLE  I 

SAMPLE  COMPARISON  OF  RATES  AND  ESTIMATES 


Pl  =  0.5,p2  = 

:0.1,T=l,d  = 

5 

Number 

Co-  Ob-  True  Raw  No  Cross- 

Cross-validation 

variate  served 

validation 

i 

2fi  Si 

ii  Xi(raw)  Xi(GA) 

Aj(Stu.t) 

Xi(GA) 

Xi(Stu.t) 

l 

+1  152 

87.57  76.00  74.56 

76.28 

53.56 

73.68 

2 

+  1  3 

1.44  1.50  1.66 

1.62 

1.66 

1.66 

3 

+1  2 

0.47  1.00  1.17 

1.23 

1.17 

126 

4 

+1  2 

3.51  1.00  1.17 

1.23 

1.17 

126 

5 

+1  0 

0.65  0.00  0.19 

0.47 

0.21 

0.52 

6 

+1  0 

0.16  0.00  0.19 

0.47 

0.21 

052 

7 

+1  5 

2.07  2.50  2.64 

2.45 

2.63 

2.48 

8 

+1  3 

2.07  1.50  1.66 

1.62 

1.66 

1.66 

9 

+1  8 

1.20  4.00  4.10 

3.76 

4.10 

3.77 

10 

+1  1 

1.11  0.50  0.68 

0.85 

0.68 

0.86 

11 

-1  9 

2.87  4.50  4.31 

4.15 

4.28 

4.16 

12 

-1  2 

0.80  1.00  1.10 

1.19 

1.10 

1.19 

13 

-1  3 

0.43  1.50  1.56 

1.59 

1.56 

1.57 

14 

-1  0 

1.22  0.00  0.18 

0.46 

0.19 

0.52 

15 

-1  9 

5.44  4.50  4.31 

4.15 

4.28 

4.16 

16 

-1  7 

3.18  3.50  3.40 

3.26 

3.38 

327 

17 

-1  11 

5.40  5.50  5.23 

5.10 

5.16 

5.08 

18 

-1  0 

0.40  0.00  0.18 

0.46 

0.19 

0.52 

19 

-1  3 

0.66  1.50  1.56 

1.59 

1.56 

1.57 

20 

-1  0 

0.03  0.00  0.18 

0.46 

0.19 

0.52 

MSE's:  7.83  9.16 

734 

58.94 

10.61 

TABLE  H 

SAMPLE  COMPARISON  OF  RATES  AND  ESTIMATES 

Pi  =  0.5,  P2  = 

0.1,  T  =  1,  d  = 

5 

Observed  True 

Raw 

No  Cross- 
validation 

Cross-validation 

Si 

A>i 

Xj(raw) 

Xj(GA) 

Xj(Stu.t) 

Xi(GA) 

Xi(Stu.t) 

267 

138.43 

133.50 

132.11 

137.56 

101.48 

130.66 

0 

0.01 

0.00 

0.17 

0.48 

0.19 

0.51 

3 

1.35 

1.50 

1.65 

1.64 

1.65 

1.65 

12 

6.55 

6.00 

6.10 

5.62 

6.10 

5.59 

0 

0.29 

0.00 

0.17 

0.48 

0.19 

0.51 

2 

2.21 

1.00 

1.16 

1.25 

1.16 

1.25 

2 

1.32 

1.00 

1.16 

1.25 

1.16 

1.25 

5 

3.47 

2.50 

2.64 

2.46 

2.64 

2.48 

0 

0.57 

0.00 

0.17 

0.48 

0.19 

0.51 

2 

0.98 

1.00 

1.16 

1.25 

1.16 

123 

7 

2.18 

3.50 

3.52 

3.42 

3.52 

339 

3 

0.47 

1.50 

1.60 

1.71 

1.60 

1.71 

36 

20.25 

18.00 

17.41 

17.48 

16.91 

17.40 

4 

1.68 

1.00 

1.12 

131 

1.12 

132 

2 

3.72 

2.00 

2.08 

2.12 

2.06 

2.12 

0 

0.39 

0.00 

0.17 

0.51 

0.18 

0.58 

6 

2.28 

3.00 

3.04 

3.97 

3.04 

2.96 

2 

1.48 

1.00 

1.12 

1.31 

1.12 

132 

9 

2.40 

4.50 

4.58 

430 

4.48 

4.29 

10 

4.46 

5.00 

4.% 

4.76 

4.95 

4.77 

MSE's 

2.22 

2.90 

1.08 

69.51 

4.09 

NOTE:  this  is  an  independent  experiment  from  the 

same  setup  as  that  of  Table  1. 
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(Row  role) 


Table  II  provides  another  illustration  of  the  estimates'  performance,  this  time 
with  fewer  extreme  outliers.  Figure  1  exhibits  histographs  of  the  mean- 
squared  error  difference,  MSE(ROLVPR)-MSE(GALVPR),  for  100  replications 
of  the  above  specific  simulation.  Note  that  the  advantage  is  in  favor  of 
ROLVPR  in  the  majority  of  the  experiments,  with  an  exceptional  advantage 
displayed  in  some  cases.  Often  it  is  in  a  few  cases  of  exceptionally  large  rates, 
and  counts,  that  ROLVPR  excels.  Figure  2  compares  the  mean-squared  errors 
of  each  shrunken  no-cross  validated  estimator  with  the  corresponding  mean 
square  error  of  the  raw-rate  estimators;  the  raw  rate  estimator  is  simply  the 
count  divided  by  tj  =  2.  For  these  data  sets  the  indications  are  that  ROLVPR 
improves  upon  RAW  most  of  the  time  when  ^2=01  and  when  {52=0.3 
(although  less  decisively),  while  RAW  improves  upon  GALVPR  most  of  the 
time;  neither  victory  is  decisive.  These  results  are  perhaps  not  surprising 
when  one  refers  to  Morris  (1983),  Theorem  1  and  subsequent  discussion.  It 
appears  that  the  convenient  conjugate  can  adapt  to  non-gamma  data  quite 
well  in  many  of  the  present  cases  at  least. 

An  undoubted  disadvantage  of  the  ROLVPR  procedure  is  its  computer 
intensivity:  computation  of  its  estimators  requires  far  more  time  than  does 
GALVPR  because  a  root  must  be  found,  (6.6),  and  a  numerical  integration 
performed.  Search  is  on  for  a  more  tractable  representation  of  a  "robust  g" 
that  permits  analytical  rather  than  computational  evaluation.  The  inverse 
Gaussian  is  a  candidate;  see  Dean  et  al.,  (1989).  Conceivably  such  an  adoption 
will  result  in  better  results  for  small-rate  situations.  Needless  to  say  RAW, 
which  quotes  Xj(RAW)  =  Sj/t;  is  by  far  the  most  economical.  Of  course  it  may 
not  be  used  if  the  covariate  value,  Xj,  changes. 


REFERENCES 

Abramowitz,  M.  and  Stegun,  I.  A.,  (1964),  Handbook  of  Mathematical 
Functions,  U.  S.  Department  of  Commerce,  National  Bureau  of 
Standards,  Applied  Mathematic  Series  55,  Seventh  Printing,  1968. 

Anscombe,  F.  J.  (1950).  Sampling  theory  of  the  negative  binomial  and 
logarithmic  series  distributions.  Biometrika  37  358-382. 

Berger,  J.  O.,  (1985),  Statistical  Decision  Theory  and  Bayesian  Analysis,  Second 
Edition,  New  York:  Springer-Verlag. 

Cox,  D.  R.,  and  Hinkley,  D.  V.,  (1974),  Theoretical  Statistics,  London: 
Chapman  and  Hall. 


18 


Dean,  C.,  Lawless,  J.  F.  and  Willmot,  G.  E.  (1989),  "A  mixed  Poisson-inverse- 
Gaussian  regression  model,"  The  Canadian  Journal  of  Statistics,  17,  171- 
181. 

Deely,  J.  J.,  and  Lindley,  D.  V.,  (1981),  "Bayes  empirical  Bayes,"  Journal  of  the 
American  Statistical  Association,  76,  833-841. 

Everitt,  B.  S.,  (1984),  An  Introduction  to  Latent  Variable  Models,  London: 
Chapman  and  Hall. 

Gaver,  D.  P.,  and  O'Muircheartaigh,  I.  G.,  (1987),  "Robust  empirical  Bayes 
analysis  of  event  rates,"  Technometrics  29, 1-15. 

Gelfand,  A.  F ,  and  Smith,  A.  F.  M.,  (1988),  "Sampling  based  approaches  to 
calculating  marginal  densities,"  Technical  Report,  Department  of 
Mathematics,  Nottingham  Statistics  Group,  University  of  Nottingham. 

Lawless,  J.  F.,  (1987a),  "Regression  methods  for  Poisson  process  data,"  Journal 
of  the  American  Statistical  Association,  82,  pp.  808-815. 

Lawless,  J.  F.,  (1987b),  Negative  binomial  regression  models.  Canad.  J.  Statist. 
15,209-226. 

McCullagh,  P.,  and  Nelder,  J.  A.,  (1983),  Generalized  Linear  Models,  London: 
Chapman  and  Hall. 

Morris,  C.  (1983),  "Parametric  empirical  Bayes  inference:  theory  and 
applications,"  (with  discussion).  Journal  of  the  American  Statistical 
Association,  78,  47-65. 

Naylor,  J.  C.,  and  Smith,  A.  F.  M.,  (1982),  "Applications  of  a  method  for  the 
efficient  computation  of  posterior  distributions,"  Applied  Statistics  31, 
214-225. 


19 


INITIAL  DISTRIBUTION  LIST 


1.  Library  (Code  0142) . 2 

Naval  Postgraduate  School 

Monterey,  CA  93943-5000 

2.  Defense  Technical  Information  Center . 2 

Cameron  Station 

Alexandria,  VA  22314 

3.  Office  of  Research  Administration  (Code  012) . 1 

Naval  Postgraduate  School 

Monterey,  CA  93943-5000 

4.  Prof.  Peter  Purdue . 1 

Code  OR- Pd 

Naval  Postgraduate  School 
Monterey,  CA  93943-5000 

5.  Department  of  Operations  Research  (Code  55) . 1 

Naval  Postgraduate  School 

Monterey,  CA  93943-5000 

6.  Prof.  Donald  Gaver,  Code  OR-Gv . ,r . 15 

Naval  Postgraduate  School 

Monterey,  CA  93943-5000 

7.  Prof.  Patricia  Jacobs . 15 

Code  OR/Jc 

Naval  Postgraduate  School 
Monterey,  CA  93943-5000 

8.  Center  for  Naval  Analyses . 1 

4401  Ford  Avenue 

Alexandria,  VA  22302-0268 

9.  Dr.  David  Brillinger . 1 

Statistics  Department 

University  of  California 
Berkeley,  CA  94720 


s 


10.  Dr.  R.  Gnanadesikan . 1 

Bellcore 

435  South  Street 
Morris  Township  NJ  07960 

11.  Prof.  Bernard  Harris . 1 

Dept,  of  Statistics 

University  of  Wisconsin 
610  Walnut  Street 
Madison,  WI  53706 

12.  Prof.  W.  M.  Hinich . 1 

University  of  Texas 

Austin,  TX  78712 

13.  Prof.  I.  R.  Savage . 1 

Dept,  of  Statistics 

Yale  University 
New  Haven,  CT  06520 

14.  Prof.  W.  R.  Schucany . 1 

Dept,  of  Statistics 

Southern  Methodist  University 
Dallas,  TX  75222 

15.  Prof.  D.  C.  Siegmund . 1 

Dept,  of  Statistics 

Sequoia  Hall 
Stanford  University 
Stanford,  CA  94305 

16.  Prof.  H.  Solomon . 1 

Department  of  Statistics 

Sequoia  Hall 
Stanford  University 
Stanford,  CA  94305 

17.  Dr.  Ed  Wegman . 1 

George  Mason  University 

Fairfax,  VA  22030 

18.  Dr.  P.  Welch . 1 

IBM  Research  Laboratory 

Yorktown  Heights,  NY  10598 


19.  Dr.  Neil  Gerr . 1 

Office  of  Naval  Research 

Arlington,  VA  22217 

20.  Prof.  Roy  Welsch . 1 

Sloan  School 

M.I.T. 

Cambridge,  MA  02139 

21.  Dr.  J.  Abrahams . 1 

Code  1111,  Room  607 


Mathematical  Sciences  Division 
Office  of  Naval  Research 
800  North  Quincy  Street 
Arlington,  VA  22217-5000 

22.  Prof.  J.  R.  Thompson . 1 

Dept,  of  Mathematical  Science 

Rice  University 
Houston,  TX  77001 

23.  Dr.  P.  Heidelberger . 1 

IBM  Research  Laboratory 

Yorktown  Heights 
New  York,  NY  10598 


24.  Prof.  M.  Leadbetter . . . 1 

Department  of  Statistics 

University  of  North  Carolina 
Chapel  Hill,  NC  27514 

25.  Prof.  D.  L.  Iglehart . 1 

Dept,  of  Operations  Research 

Stanford  University 
Stanford,  CA  94350 

26.  Prof.  J.  B.  Kadane . 1 

Dept,  of  Statistics 

Carnegie-Mellon  University 
Pittsburgh,  PA  15213 


$ 


27.  Prof.  J.  Lehoczky . 1 

Department  of  Statistics 

Carnegie-Mellon  University 
Pittsburgh,  PA  15213 

28.  Dr.  J.  Maar  (R513) . 1 

National  Security  Agency 

Fort  Meade,  MD  20755-6000 

29.  Prof.  M.  Mazumdar . 1 

Dept,  of  Industrial  Engineering 

University  of  Pittsburgh 
Pittsburgh,  PA  15235 

30.  Prof.  M.  Rosenblatt . 1 

Department  of  Mathematics 

University  of  California,  San  Diego 
La  Jolla,  CA  92093 

31.  Prof.  H.  Chernoff . : . 1 

Department  of  Statistics 

Harvard  University 
1  Oxford  Street 
Cambridge,  MA  02138 

32.  Dr.T.J.  Ott . 1 

Bellcore 

435  South  Street 

Morris  Township,  NJ  07960 

33.  Dr.  Alan  Weiss . 1 

AT&T  Bell  Laboratories 

Mountain  Avenue 
Murray  Hill,  NJ  07974 

34.  Prof.  Joseph  R.  Gani . 1 

Mathematics  Department 

University  of  California 
Santa  Barbara,  CA  93106 

35.  Prof.  Frank  Samaniego . 1 

Statistics  Department 

University  of  California 
Davis,  CA  95616 


36. 


1 


Dr.  James  McKenna . 

AT&T  Bell  Laboratories 
Mountain  Avenue 
Murray  Hill,  NJ  07974 

37.  Commander . 1 

Attn:  G.  F.  Kramer 

Naval  Ocean  Systems  Center 
Code  421 

San  Diego,  CA  92152-5000 

38.  Prof.  Tom  A.  Louis . 1 

School  of  Public  Health 

University  of  Minnesota 
Mayo  Bldg.  A460 
Minneapolis,  MN  55455 

39.  Dr.  Nan  Laird . 1 

Biostatistics  Dept. 

Harvard  School  of  Publfc  Health 
677  Huntington  Ave. 

Boston,  MA  02115 

40.  Dr.  Marvin  Zelen . 1 

Biostatistics  Department 

Harvard  School  of  Public  Health 
677  Huntington  Ave. 

Boston,  MA  02115 

41.  Dr.  John  Ora v . 1 

Biostatistics  Department 

Harvard  School  of  Public  Health 
677  Huntington  Ave. 

Boston,  MA  02115 

42.  Prof.  R.  Douglas  Martin . 1 

Department  of  Statistics,  GN-22 

University  of  Washington 
Seattle,  WA  98195 
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43.  Prof.  W.  Stuetzle . 

Department  of  Statistics 
University  of  Washington 
Seattle,  WA  98195 

44.  Prof.  F.  W.  Mosteller . 1 

Department  of  Statistics 

Harvard  University 
1  Oxford  St. 

Cambridge,  MA  02138 

45.  Dr.  D.  C.  Hoaglin . 1 

Department  of  Statistics 

Hardvard  University 
1  Oxford  Street 
Cambridge,  MA  02138 


46.  Prof.  N.  D.  Singpurwalla . 1 

George  Washington  University 

Washington,  DC  20052 

47.  Center  for  Naval  Analysis . . . 1 

2000  Beauregard  Street 

Alexandria,  VA  22311 

48.  Prof.  George  S.  Fishman . . . 1 


Curr.  in  OR  &  Systems  Analysis 
University  of  North  Carolina 
Chapel  Hill,  NC  20742 

49.  Dr.  Alan  F.  Petty . 1 

Code  7930 

Naval  Research  Laboratory 
Washington,  DC  20375 

50.  Prof.  Bradley  Efron . 1 

Statistics  Dept. 

Sequoia  Hall 
Stanford  University 
Stanford,  CA  94305 


s 


1 


51.  Prof.  Carl  N.  Morris . 

Statistics  Department 
Harvard  University 
1  Oxford  St. 

Cambridge,  MA  02138 


52.  Dr.  John  E.  Rolph . 1 

RAND  Corporation 

1700  Main  St. 

Santa  Monica,  CA  90406 

53.  Prof.  Linda  V.  Green . 1 


Graduate  School  of  Business 
Columbia  University 
New  York,  NY  10027 

54.  Dr.  David  Burman . 1 

AT&T  Bell  Telephone  Laboratories 

Mountain  Avenue 
Murray  Hill,  NJ  07974  ' 

55.  Dr.  Ed  Coffman . 1 

AT&T  Bell  Telephone  Laboratories 

Mountain  Avenue 
Murray  Hill,  NJ  07974 

56.  Prof.  William  Jewell . 1 

Operations  Research  Department 

University  of  California,  Berkeley 
Berkeley,  CA  94720 

57.  Prof.  D.  C.  Siegmund . 1 

Dept,  of  Statistics 

Sequoia  Hall 
Stanford  University 
Stanford,  CA  94305 

58.  Operations  Research  Center,  Rm  E40-164 . 1 

Massachusetts  Institute  of  Technology 

Attn:  R.  C.  Larson  and  J.  F.  Shapiro 
Cambridge,  MA  02139 
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59.  Arthur  P.  Hurter,  Jr . 1 

Professor  and  Chairman 

Dept,  of  Industrial  Engineering  and  Management  Sciences 
Northwestern  University 
Evanston,  IL  60201-9990 

60.  Institute  for  Defense  Analysis . 1 

1800  North  Beauregard 

Alexandria,  VA  22311 

61.  Prof.  J.  W.  Tukey . 1 

Statistics  Dept.,  Fine  Hall 

Princeton  University 
Princeton,  NJ  08540 

62.  Dr.  Daniel  H.  Wagner . 1 

Station  Square  One 

Paoli,  PA  19301 

63.  Dr.  Colin  Mallows . . . 1 

AT&T  Bell  Telephone  Laboratories 

Mountain  Avenue 
Murray  Hill,  NJ  07974 

64.  Dr.  D.  Pregibon . 1 

AT&T  Bell  Telephone  Laboratories 

Mountain  Avenue 

Murray  Hill,  NJ  07974  '' 

65.  Dr.  Jon  Kettenring . . . 1 

Bellcore 

435  South  Street 

Morris  Township,  NJ  07960 

66.  Prof.  David  L.  Wallace . 1 

Statistics  Dept. 

University  of  Chicago 
5734  S.  University  Ave. 

Chicago,  IL  60637 

67.  Dr.  S.  R.  Dalai . 1 

AT&T  Bell  Telephone  Laboratories 

Mountain  Avenue 
Murray  Hill,  NJ  07974 


1 


68.  Dr.  M.  J.  Fischer . 

Defense  Communications  Agency 
1860  Wiehle  Avenue 
Reston,  VA  22070 


69.  Dr.  Prabha  Kumar 

Defense  Communications  Agency . 1 

1860  Wiehle  Avenue 
Reston,  VA  22070 

70.  Dr.  B.  Doshi . 1 

AT&T  Bell  Laboratories 

HO  3M-335 
Holmdel,  NJ  07733 

71.  Dr.  D.  M.  Lucantoni . 1 

AT&T  Bell  Laboratories 

Holdmel,  NJ  07733 

72.  Dr.  V.  Ramaswami . ;. . 1 

MRE  2Q-358 


Bell  Communications  Research,  Inc. 

435  South  Street 
Morristown,  NJ  0 7960 

73.  Prof.  G.  Shantikumar . 1 

The  Management  Science  Group 

School  of  Business  Administration 
University  of  California 
Berkeley,  CA.  94720 

74.  Dr.  D.  F.  Daley . 1 

Statistic  Dept.  (I.A.S.) 

Australian  National  University 
Canberra,  A.C.T.  2606 
AUSTRALIA 

75.  Dr.  Guy  Fayolle . 1 

I.N.R.I.A. 

Dom  de  Voluceau-Rocquencourt 
78150  Le  Chesnay  Cedex 
FRANCE 


76.  Professor  H.  G.  Daellenbach . 1 

Department  of  Operations  Research 

University  of  Canterbury 
Christchurch,  NEW  ZEALAND 

77.  Koh  Peng  Kong . 1 

OA  Branch,  DSO 

Ministry  of  Defense 
Blk  29  Middlesex  Road 
SINGAPORE  1024 

78.  Professor  Sir  David  Cox . 1 

Nuffield  College 

Oxford,  OXI  INF 
ENGL  AND 

79.  Dr.  A.  J.  Lawrence . 1 

Department  of  Mathematical  Statistics 

University  of  Birmingham 
P.  O.  Box  363 
Birmingham  B15.  2TT 
ENGLAND 

80.  Dr.  F.  P.  Kelly . 

Statistics  Laboratory 
16  Mill  Lane 
Cambridge 
ENGLAND 

81.  Dr.  R.  J.  Gibbens . . . 1 

Statistics  Laboratory 

16  Mill  Lane 

Cambridge 

ENGLAND 

82.  Dr.  John  Copas . 1 

Dept,  of  Mathematical  Statistics 

University  of  Birmingham 
P.  O.  Box  363 
Birmingham  B15  2TT 
ENGLAND 
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83.  Dr.  D.  Vere-Jones . 

Dept,  of  Mathematics 

Victoria  University  of  Wellington 
P.  O.  Box  196 
Wellington 
NEW  ZEALAND 

84.  Prof.  Guy  Latouche . 

University  Libre  Bruxelles 

C.  P.  212 

Blvd.  De  Triomphe 
B-1050  Bruxelles 
BELGIUM 


